/*#if !defined(NDEBUG)
  #define BOOST_MULTI_INDEX_ENABLE_INVARIANT_CHECKING
  #define BOOST_MULTI_INDEX_ENABLE_SAFE_MODE
  #endif
*/
#ifdef _WIN32
#include <conio.h>
#endif



#include <ctype.h>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <iterator>
#include <string>
#include <time.h>

#include <set>
#include <cmath>
#include <algorithm>
#include <map>
#include <vector>
#include <iostream>
#include <iterator>
#include <functional>
#include <algorithm>
#include <valarray>
#include <boost/foreach.hpp>
#include <boost/timer.hpp>
#include <boost/progress.hpp>
#include <boost/multi_index_container.hpp>
#include <boost/multi_index/member.hpp>
#include <boost/multi_index/ordered_index.hpp>
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics.hpp>
#include <boost/accumulators/statistics/tail_quantile.hpp>




#ifdef USE_ARCHIVE
#include <boost/archive/tmpdir.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/serialization/base_object.hpp>
#include <boost/serialization/utility.hpp>
#include <boost/serialization/list.hpp>
#include <boost/serialization/assume_abstract.hpp>
#endif




#include "GenThreads.h"

#include "loader.h"
#include "IdTest.h"


using boost::multi_index_container;
using namespace boost::multi_index;

///////////////////////////////
#define stop exit(0);
///////////////////////////////
  /* tags for accessing the corresponding indices of particle_set */

  struct id{};
  //	struct idu{};
  struct snap{};
  //	struct snapu{};
#ifdef R_INDEX
  struct R{};
#endif
////////////////////////////////
class CAnalysis:public CGenThreads{
 public:
  CAnalysis();
  ~CAnalysis(){delete m_ComVec;delete m_ApVec;};

  /* an particle record holds its ID, name and age */
struct particle
  {
    friend class boost::serialization::access;
    // When the class Archive corresponds to an output archive, the
    // & operator is defined similar to <<.  Likewise, when the class Archive
    // is a type of input archive the & operator is defined similar to >>.
    template<class Archive>
    void serialize(Archive & ar, const unsigned int version)
    {
      ar & id;
      ar & snap;
      ar & Pos[0];
      ar & Pos[1];
      ar & Pos[2];
      ar & Vel[0];
      ar & Vel[1];
      ar & Vel[2];
      ar & Rho;
      ar & Hsml;
#ifdef R_INDEX
      ar & R;
      //			 ar & RR[0];
      //			 ar & RR[1];
#endif
    }
    int         id;
    unsigned short snap;
    float Pos[3];
    float Vel[3];
    float Rho;
    float Hsml;
#ifdef R_INDEX
    float R;
#endif
    particle(int id_):id(id_){
      Pos[0]=0;
      Pos[1]=0;
      Pos[2]=0;
      Vel[0]=0;
      Vel[1]=0;
      Vel[2]=0;
      Rho=0;
      Hsml=0;
#ifdef R_INDEX
      R=0;
#endif	
    };
    particle(){};
    particle(int id_,unsigned short snap_,float *Pos_, float *Vel_, 
	     float _Rho=0, float _Hsml=0):
      id(id_),snap(snap_),Rho(_Rho),Hsml(_Hsml)
    {
      Pos[0]=Pos_[0];
      Pos[1]=Pos_[1];
      Pos[2]=Pos_[2];
      Vel[0]=Vel_[0];
      Vel[1]=Vel_[1];
      Vel[2]=Vel_[2];		
#ifdef R_INDEX
      R= GetR();
#endif
    };
    public:
    particle& operator*=(const double& right)
    {
      this->Pos[0] *= (float)right;
      this->Pos[1] *= (float)right;
      this->Pos[2] *= (float)right;
      this->Vel[0] *= (float)right;
      this->Vel[1] *= (float)right;
      this->Vel[2] *= (float)right;
      this->Hsml *= (float)right;
      this->Rho *= (float)right;
      return *this;
    }
    /////////////////////////////////
    //make a mean position for a two particles
    particle& operator+=(const particle& right)
    {
      this->Pos[0] += right.Pos[0];
      this->Pos[1] += right.Pos[1];
      this->Pos[2] += right.Pos[2];
      this->Vel[0] += right.Vel[0];
      this->Vel[1] += right.Vel[1];
      this->Vel[2] += right.Vel[2];
      this->Hsml += right.Hsml;
      this->Rho += right.Rho;          
    return *this;
    }
    /////////////////////////////////
    float  GetR(void) const{return std::sqrt(Pos[0]*Pos[0]+Pos[1]*Pos[1]+Pos[2]*Pos[2]);};

    friend std::ostream& operator<<(std::ostream& os, const particle& e)
    {
      //os<<"SNAP: "<<e.snap<<" ID: "<<e.id<<" R: "<<e.R<<std::endl;
      os<<e.snap;
      os.width(16);
      os<<"\t"<<e.id;
      os.precision(4);
      os.width(16);
      os.setf( std::ios::fixed, std::ios::floatfield ) ;
      os<<" "<<e.Pos[0];
      os.width(16);
      os<<"\t"<<e.Pos[1];
      os.width(16);
      os<<"\t"<<e.Pos[2];
      os.width(10);			
      os<<"\t"<<e.GetR()<<std::endl;

      return os;
    }
    bool operator<(const particle& e)const{return GetR()<e.GetR();}

    void  get_velz(float *A)const//Ap, Ar, Az
    {
      float x=Pos[0];
      float y=Pos[1];
      float z=Pos[2];
      float Ax=Vel[0];
      float Ay=Vel[1];
      A[2]=Vel[2];

      float R=std::sqrt(x*x+y*y);
      float rr=GetR();
      
      float cos_phi=x/R, sin_phi=y/R,
	cos_theta=z/rr, sin_theta=R/rr;

      A[1] = Ax*cos_phi +Ay*sin_phi;
      A[0] = -Ax*sin_phi+Ay*cos_phi;

    }

  };
  struct comp_id
  {
    // compare an ID and a particle
    bool operator()(int x,const particle& e2)const{return x<e2.id;}

    // compare an employee and an ID
    bool operator()(const particle& e1,int x)const{return e1.id<x;}
  };
  /* see Compiler specifics: Use of member_offset for info on
   * BOOST_MULTI_INDEX_MEMBER
   */

  /* Define a multi_index_container of particles with following indices:
   *   - a unique index sorted by particle::int,
   *   - a non-unique index sorted by particle::name,
   *   - a non-unique index sorted by particle::R.
   */
  /* change ordered_non_unique to ordered_unique for ID to avoid a partcile with
   * the same ID
   */

  typedef multi_index_container<
    particle,
    indexed_by<
    ordered_non_unique<
    tag<id>,  BOOST_MULTI_INDEX_MEMBER(particle,int,id)>,
    ordered_non_unique<
    tag<snap>,BOOST_MULTI_INDEX_MEMBER(particle,unsigned short,snap)>
#ifdef R_INDEX
    ,ordered_non_unique<
	  
    tag<R>, BOOST_MULTI_INDEX_MEMBER(particle,float,R)>
#endif
    >
	  
    > particle_set;


  typedef multi_index_container<
    int,
    indexed_by<
    ordered_unique<identity<int> >
    >
    > id_unique_set;

  typedef multi_index_container<
    string,
    indexed_by<
    ordered_unique<identity<string> >
    >
    > snap_unique_set;
  typedef multi_index_container<
    double,
    indexed_by<
    ordered_unique<identity<double> >
    >
    > time_unique_set;

int  parallel_hello_data(unsigned short  isnap, unsigned short myThreadID);	
  template<typename Tag,typename MultiIndexContainer>
    void print_out_by(
		      const MultiIndexContainer& s,
		      Tag* =0 /* fixes a MSVC++ 6.0 bug with implicit template function parms */
		      )
    {
      /* obtain a reference to the index tagged by Tag */

      const typename boost::multi_index::index<MultiIndexContainer,Tag>::type& i=
	get<Tag>(s);

      typedef typename MultiIndexContainer::value_type value_type;

      /* dump the elements of the index to cout */

      std::copy(i.begin(),i.end(),std::ostream_iterator<value_type>(std::cout));
    }

	
  ////////////////////////////////////////////////////////////////////////
#ifdef DEBUG_ASS
  void print_out_by_id()
    { print_out_by<id>(ps);}
  void print_out_by_snap()
    { print_out_by<snap>(ps);}
#endif
  ////////////////////////////////////////////////////////////////////////
  int load_data(string filename, int TYPE,string snapio);
  void SetupCom(int a, int b);
  vector<float> FillAp(int isnap);
  int parallel_load_data(unsigned short isnap, bool append=false,unsigned short myThreadID=0);
  void DumpOneOrbit( int inID);
  void OpenOrbitFile(int myID);
  void DumpOrbits();
  bool Load(void);
  void Save(void);
  CComPos GetCom(int isnap);  
  bool isInSel(int inID);
#ifndef DEBUG_ASS
  int GetNTracked(){return idu_sel.size();};
  int GetNPs(){return ps.size();};
#endif
  void DumpOrbitsOAF();
  //////////////////////////////////////////////////////////////////////
  void my_print_stats();
 protected:
  // create and open a character archive for output
  std::ofstream ofs;

  //vector<int> m_type;/it is unclear why we need it?

  particle_set ps;
  id_unique_set idu;
  snap_unique_set snapu;
  time_unique_set timeu;
  /************************/
  id_unique_set idu_sel;
  /********************************/
  particle_set::nth_index<0>::type& sort_by_id;
  particle_set::nth_index<1>::type& sort_by_snap;
#ifdef R_INDEX
  particle_set::nth_index<2>::type& sort_by_R;
#endif
  /********************************/
  void WriteMeanFile(vector<particle> data,string meanfilename);
  bool InitData(int isnap, CGetData *&data);
  //functions doing real analysis job
  void parallel_profiler_data(int isnap);
  void parallel_Test_data(int isnap);
//////////////////////////////////////////////
  public:
  //////////////////////////////////////////////////
  void EmptyIDsContainers(){
    ps.clear();
    idu.clear();
    snapu.clear();
    timeu.clear();
    // m_type.clear();
};
  //////////////////////////////////////////////////
  /// helper functions for real job functions;
	  void do_parallel_load_on_SNAPS(int a, int b);
	  void do_parallel_Profiler_on_SNAP(int a, int b);
	  void do_parallel_Test_on_SNAP(int a, int b);	 
	  void do_parallel_save_on_OAF(int a, int b);
//////////////////////////////////////////////
 protected:
  int count;
  CComVec *m_ComVec;// this keeps COM per snapshot
  CApVec *m_ApVec;//this keeps Pericenter,Apocenter peaks per snapshot
  boost::mutex m_ps_fill_mutex;
  ////////////////////////////////
  std::vector<unsigned int> sindex;// Used for sorting the data by Radius;
  ///////////////////////////
  //Ap stuff
  vector<double> xAp, Ap;
  std::pair<double, double>   BarRot;
  string m_OmegaFile;
  string 	fnameApPeaks;
  CIdTest id_set;
};


